home *** CD-ROM | disk | FTP | other *** search
/ Mac Magazin/MacEasy 21 / Mac Magazin and MacEasy Magazine CD - Issue 21.iso / Wissenschaft & Technik / yorick12vr1-nofpu folder / include / series.i < prev    next >
Text File  |  1995-07-26  |  5KB  |  163 lines
  1. /*
  2.    SERIES.I
  3.    Routines for handling geometric series (e.g.- tapered zoning).
  4.  
  5.    $Id: series.i,v 1.1 1993/08/27 18:50:06 munro Exp $
  6.  */
  7. /*    Copyright (c) 1994.  The Regents of the University of California.
  8.                     All rights reserved.  */
  9.  
  10. func series_s(r, n)
  11. /* DOCUMENT series_s(r, n)
  12.      returns the sum s of the finite geometric series
  13.         1 + r + r^2 + r^3 + ... + r^n
  14.      Using n<0 is equivalent to using the reciprocal of r, that is,
  15.         series_s(r, -n) == series_s(1./r, n)
  16.      R or N or both may be arrays, as long as they are conformable.
  17.    SEE ALSO: series_r, series_n
  18.  */
  19. {
  20.   /* force conformability immediately */
  21.   dims= dimsof(r, n);
  22.   rr= r + array(0.0, dims);
  23.   nn= n + array(0, dims);
  24.  
  25.   /* form result array, initialized to n==0 result */
  26.   s= array(1.0, dims);
  27.  
  28.   /* subdivide into n==0, n>0, n<0, and r==1.0 cases */
  29.   rnot= rr!=1.0;
  30.  
  31.   mask= (nn>0 & rnot);   /* n>0 case */
  32.   list= where(mask);
  33.   if (numberof(list)) {
  34.     rrx= rr(list);
  35.     s= merge((rrx^(1+nn(list))-1.0)/(rrx-1.0), s(where(!mask)), mask);
  36.   }
  37.  
  38.   mask= (nn<0 & rnot);   /* n<0 case */
  39.   list= where(mask);
  40.   if (numberof(list)) {
  41.     rrx= 1.0/rr(list);
  42.     s= merge((rrx^(1-nn(list))-1.0)/(rrx-1.0), s(where(!mask)), mask);
  43.   }
  44.  
  45.   list= where(!rnot);   /* r==1.0 case */
  46.   if (numberof(list)) s= merge(s(where(rnot)), abs(nn(list)), rnot);
  47.  
  48.   return s;
  49. }
  50.  
  51. func series_r(s, n)
  52. /* DOCUMENT series_r(s, n)
  53.      returns the ratio r of the finite geometric series, given the sum s:
  54.         1 + r + r^2 + r^3 + ... + r^n = s
  55.      Using n<0 will return the the reciprocal of n>0 result, that is,
  56.         series_r(s, -n) == 1.0/series_r(s, n)
  57.      If n==0, returns s-1 (the n==1 result).
  58.      S or N or both may be arrays, as long as they are conformable.
  59.    SEE ALSO: series_s, series_n
  60.  */
  61. {
  62.   /* force conformability immediately */
  63.   dims= dimsof(s, n);
  64.   ss= s + array(0.0, dims);
  65.   nn= n + array(0, dims);
  66.  
  67.   /* form result array, initialized to abs(n)<2 result */
  68.   r= ss-1.0;
  69.  
  70.   /* work only with n>0 -- take reciprocals at end if necessary */
  71.   nneg= nn<0;
  72.   nn= abs(nn)+1;
  73.   nbig= nn>2;
  74.  
  75.   /* compute an approximate result which has exact values and
  76.      derivatives at s==1, s==n, and s->infinity --
  77.      different approximations apply for s>n and s<n */
  78.   mask= nbig & (ss>nn);
  79.   list= where(mask);
  80.   if (numberof(list)) {
  81.     sx= ss(list);
  82.     nx= nn(list);
  83.     pow= 1.0/(nx-1.0);
  84.     npow= nx^pow - 1.0;
  85.     n2r= 1.0/(nx-2.0);
  86.     A= (2.0-nx*npow)*n2r;
  87.     B= (2.0*npow-nx*pow)*nx*n2r;
  88.     r= merge(sx^pow - pow + A*(nx/sx)^pow + B/sx, r(where(!mask)), mask);
  89.   }
  90.   mask= nbig & (ss<=nn);
  91.   list= where(mask);
  92.   if (numberof(list)) {
  93.     sx= ss(list);
  94.     nx= nn(list);
  95.     sn= (sx-1.0)/(nx-1.0);
  96.     n2r= 1.0/(nx*nx);
  97.     r= merge(1.0 - 1.0/sx + n2r*sn*sn*(nx+1.0 - sn), r(where(!mask)), mask);
  98.   }
  99.  
  100.   /* Polish the approximation using Newton-Raphson iterations.
  101.      There are never many of these, so do the entire vector together.  */
  102.   mask= nbig & (ss!=nn);
  103.   list= where(mask);
  104.   if (numberof(list)) {
  105.     rx= r(list);
  106.     ss= ss(list);
  107.     nn= nn(list);
  108.     for (;;) {
  109.       rr= rx-1.0;
  110.       rn= rx^(nn-1);
  111.       rrss= rr*ss;
  112.       delta= rrss - (rx*rn-1.0);
  113.       if (allof(abs(delta)<=1.e-9*abs(rrss))) break;
  114.       rx+= delta/(nn*rn-ss);
  115.     }
  116.     /* try to get it to machine precision */
  117.     if (anyof(delta)) rx+= delta/(nn*rn-ss);
  118.     r= merge(rx, r(where(!mask)), mask);
  119.   }
  120.  
  121.   list= where(nneg);
  122.   if (numberof(list)) r= merge(1.0/r(list), r(where(!nneg)), nneg);
  123.  
  124.   return r;
  125. }
  126.  
  127. func series_n(r, s)
  128. /* DOCUMENT series_n(r, s)
  129.      returns the minimum number n of terms required for the geometric
  130.      series
  131.         1 + r + r^2 + r^3 + ... + r^n = s
  132.      to reach at least the given value s.  An alternate viewpoint is
  133.      that n is the minimum number of terms required to achieve the
  134.      sum s, with a ratio no larger than r.
  135.      Returns 0 if r<1 and s>1/(1-r), or if s<1.
  136.      The routine makes the most sense for r>1 and s substantially
  137.      greater than 1.  The intended use is to determine the minimum
  138.      number of zones required to span a given thickness t with a given
  139.      minimum zone size z, and maximum taper ratio r (assumed >1 here):
  140.         n= series_n(r, t/z);
  141.      With this n, you have the option of adjusting r or z downwards
  142.      (using series_r or series_s, respectively) to achieve the final
  143.      desired zoning.
  144.      R or S or both may be arrays, as long as they are conformable.
  145.    SEE ALSO: series_s, series_r
  146.  */
  147. {
  148.   n= 1.0 + (r-1.0)*s;
  149.   bad= (n<=0.0 | s<1.0);
  150.   list= where(bad);
  151.   if (numberof(list))
  152.     n= merge(array(1.0, numberof(list)), n(where(!bad)), bad);
  153.   mask= r==1.0 & !bad;
  154.   list= where(mask);
  155.   if (numberof(list))
  156.     n= merge((s+0.0*n)(list), n(where(!mask)), mask);
  157.   mask= !mask & !bad;
  158.   list= where(mask);
  159.   if (numberof(list))
  160.     n= merge(log(n(list))/log(r(list)), n(where(!mask)), mask);
  161.   return long(ceil(n))-1;
  162. }
  163.